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Abstract 

Qj We analyze the effect of increasing charge density on the Fixed Node Errors in 

■^ Diffusion Monte Carlo by comparing FN-DMC calculations of the total ground state 

G energy on a 4 electron system done with a Hartree-Fock based trial wave function 

(— I to calculations by the same method on the same system using a Configuration 

O Interaction based trial wave function. We do this for several different values of 

c/3 nuclear charge, Z. The Fixed Node Error of a Hartree-Fock trial wave function for 

• ^H a 4 electron system increases linearly with increasing nuclear charge. 
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1 Introduction 



1.1 place of QMC amongst other electronic structure models 



Quantum Monte Carlo (QMC) represents a promising methodology for solv- 

.^ ing electronic structure and quantum many-body problems in general. The 

/\^ common thread of QMC methods is the use of stochastic sampling both to 

5^ evaluate expectation values for many-body wave functions as well as to solve 

the stationary Schrodinger equation [Ij. In addition to being a useful tool for 

treating interactions in a many-body framework, QMC methods can handle 

many different systems ranging from atoms to molecules to extended crystals. 

The obtained results for energy differences are typically within a few percent 

of experimental values. Even on absolute scale QMC correlation energies are 
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remarkably accurate, reaching 90-95% of the exact values. The favorable scal- 
ing in the number of particles makes QMC very promising both in applications 
and also in further theoretical developments. 

The key approximation which hampers further increase in accuracy is the 
so-called fixed-node approximation which is used to overcome the notorious 
fermion sign problem. The fermionic antisymmetry implies that wave func- 
tions will have both positive and negative negative values; however, the sta- 
tistical approaches such as QMC typically require non-negative distributions 
otherwise they become inefficient. One means of overcoming this fundamental 
difficulty is to use the node (wave function zero locus) between the positive 
and the negative regions as a boundary condition in solving the Schrodinger 
equation. This avoids the sign problems since one finds the absolute value 
of the wave function with prescribed boundaries. In the case that one would 
know the exact nodal hypersurface, the method would provide the exact solu- 
tion. Unfortunately one seldom possesses such knowledge. Therefore we have 
to use the best available nodes which come from optimized trial (or varia- 
tional) wave functions. We fix the nodal surface of the solution to be identical 
to the nodes of the trial function and solve the Schrodinger equation with 
approximate nodes. The fixed-node approximation provides an upper bound 
to the true energy, and in practice it provides very remarkable accuracy as 
mentioned above. 

The goal of this work is to study the impact of increasing electronic density 
on the fixed-node bias. This requires a system which shows significant error 
in the mean-field (Hartree-Fock) nodes but for which we know an excellent 
approximation to the exact node [2|5B ] which is not too difficult to construct 
across the range of densities. In particular, we focus on studying a "Be-like" 
series of four electron systems in a Coulombic potential with varying nuclear 
charge Z = 3 to Z = 28. We contrast the fixed-node bias of the Hartree- 
Fock vs. two-configuration trial wave functions and analyze the corresponding 
trends as Z increases. 



1.2 Basics of DMC and the fixed-node approximation. 



Let us briefiy mention the key facts about the fixed-node diffusion Monte 
Carlo method (FN-DMC). Consider an operator exp{—TH) where H is the 
Hamiltonian and r is a real parameter (imaginary time) . It is straightforward 
to show that applying this operator to an initial trial wave function \l/r 

lim exp(-rif)\l/T = v&o 

T— S>00 

projects out the ground state \I'o within the given symmetry class. This pro- 
jection can be formulated as a solution of the imaginary-time Schrodinger 



equation 

/(R, t + t') = J G(R f- R', r')/(R, r) 

where 

G'(Rf-R',r') = ^T(R)^?^(R')(R|e""^|R') 

is the Green function and R is the vector of electron coordinates. We have 
introduced the importance samphng by \l/ so that the desired ground state $ 
is obtained from the solution /(R, r) = $(R, r)\E'7-(R) for large r. Note that 
the fixed-node condition implies that $ vanishes at the nodes of "^t and also 
that /(R, r) > for any R, r. Further details about the FN-DMC method 
can be found in Ref. p. 



1.3 Origin of nodal errors: topology of two versus four nodal domains in ^e 
system. 



Let us consider the lowest state of four electrons in the Coulomb potential 
with charge Z (we assume atomic units throughout). To aid in understanding 
the nature of fixed-node errors for this system, we recall the properties of the 
nodal domain topology of the two different trial wave functions. The Hartree- 
Fock (HF) trial wave function for Be is given by a Slater determinant which 
is block diagonal in spin so that it can be broken into a product of the spin 
channels 

^hf(R) =det[0is(ri),02s(r2)] x det[0i^(r3), 02^(^4)] 
where rj = |rj|. The nodes are defined implicitly for any trial wave function as 

^t(R) = 

Let us analyze the determinant for the non-interacting orbitals first. We can 
write analytically the hydrogenic functions as 

01. = e-^' 



(.2. = [1 - {Zr/2)]e 



-Zr/2 



and 

so that we get 

h{ri)h{r2)[g{ri) - g{r2)] = 

where h{r) is a function which does not vanish for any finite r while 

g{r) = exp(-Zr/2) - 1 + {Zr/2) 

It is then easy to show that the function g{r) is monotonous since g'{r) > 
for any r > 0. Using numerically accurate HF orbitals for Be, it is possible 



to carry out similar arrangements and to demonstrate qualitatively the same 
behavior. The monotonicity is important since it immediately implies that 
equation [g{ri) — g{r2)] = is equivalent to ri — r2 = and also that this is 
the only solution. The node of the HF/non-interacting trial function is then 
given by 

(ri -r2)(r3-r4) = 0, 

which clearly shows that there are 2x2 = 4 nodal pockets [2]. However, 
it has been found some time ago that the correct number of nodal domains 
is two, see, for example, Ref. [2]. (For our purposes it is very instructive 
that the corresponding fixed- node bias in the correlation energy is significant, 
about 10% of the correlation energy |3J.) The accurate nodal surface for the ^S 
ground state is actually remarkably well described by two configurations wave 
function where HF is augmented by adding 2s^ — )■ 2p^ double excitation which 
corresponds to a near-degeneracy effect. This configuration state function is 
a sum of three Slater determinant products det[0is, 02pJ det[0is, 02pJ for i = 
x,y,z as required by the ^S symmetry. The 2-configuration wave function is 
thus given as 

where the d^ and di are the variational coefficients. 



2 Trial wave function 

We write the many-body wave-function as a product of an antisymmetric 
determinantal part, "^a, times a Jastrow correlation factor 

Our Jastrow factor includes 1-, 2-, and 3-body terms and has the form 

ilk 

' / , ^mOm\'ij) 
ijm 

+ Yl Ckim\ak{rii)ai{rji) + ak{rji)ai{rii)>b^{rij) 

ijlklm 

where the a^ are one-body basis terms, the bm are two-body basis terms, 
the Cfc, Cm, Ckim are coefficients, and the i,j label electrons while capital / 
labels ions. The sums in the Jastrow factor can be understood as taking into 
account electron-ion, electron-electron, and electron-electron-ion correlation, 
respectively [T]. Further details about the correlation functions can be found 
in Ref.[S]. As explained above, we employ two types of trial functions. The 



first one has the HF nodes so that 
while the second state 

*A = *2con/ 

Our one-particle orbitals were generated using a numerical HF code. The 
variational parameters to be optimized in variational Monte Carlo are the 
Jastrow coefficients, Ck,Cm, and Ckim, a parameter in the basis terms a^ and bm, 
and the determinantal weight, di. These were optimized by algorithms which 
minimize the combination of energy and its variance as described elsewhere 
Ref . [5] . All QMC calculations were done using the software package QWalk [6] . 

2. 1 Dependence of Ehf o-nd Ecorr on Z 

The correlation energy is defined as customary 

E,{N, Z) = E{N, Z) - Ehf{N, Z) 

where E{N, Z) is the total non-relativistic energy of an atom with nuclear 
charge Z and A^ electrons, while Ehf{N, Z) is its Hartree-Fock counterpart. 
E{N, Z)/Z^ can be expanded in a Laurent series as in [7j 

E{N, Z) =Bo{N)Z^ + B^{N)Z + B^iN) 

+ Bs{N)Z-^ + Bi{N)Z-^ + ... 

We can use M0ller-Plesset perturbation theory to analyze the first few terms. 
In particular, one finds that E^f = Bq{N)Z^ + Bi{N)Z. The ffist term is a 
sum of the energies of occupied orbitals, 

BoiN)Z' = J2ea 

a 

where each orbital energy ea is proportional to Z^. The next term is equal to 
ffist order correction to this energy with the perturbing potential being the 
difference between the exact potential and the effective Hartree-Fock potential, 

Vpert = Vcoulomb — VhF, 

B,{N)Z = i^HFlVpertl^HF) 

which will be linear in Z. For the four electron system, Ec will increase lin- 
early with 2' [7j [8]. Chakravorty et. al. explain that this is due to the near- 
degeneracy effect since the single configuration Hartree-Fock wave function 
is not the correct zeroth order wave function in the non-interacting limit 
Z^^ = 0. This deficiency means {'^HF\^pert\'^HF) is not equal to Bi. As such. 
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Fig. 1. A comparison of the FN-DMC error for different wave functions calculated 
using values in Table [TJ The squares correspond to the HF nodes while the circles 
correspond to the 2-configuration nodes. The linear fit to the error from the HF 
nodal structure has a slope of 0.0111(1). The error bars are much smaller than the 
plot symbols. 

expansions of the HF energy and the actual total energy in Z~^ for Be and 
Be-like atoms differ in the linear coefficient. This results in a leading term 
proportional to Z in the expression for E^orr [Zj- 



3 FNDMC results and discussion 



We carried out the fixed-node calculations for trial functions with '^hf and 
^2conj nodes for Z = 3,4,5,10,12,20,28. The results are listed in Table [l] 
and plotted in Fig. [l] The linear increase with growing Z is very clear and 
is quite striking considering that it covers the range from Li~ through Be 
to very highly ionized cases. Interestingly, this linear increase matches the 
linear increase in the correlation energy based on perturbation analysis above. 
In fact, the observed slope of -0.111(1) is very close the analytically derived 
result of -0.0117 for the linear correction of the total energy in the MP2 theory 



Table 1 

FN-DMC ground state energies for ^hf and ^2conf wave functions compared to 
the estimated exact energies estimated from experiments [Tj for Z=4 through 28 
and extrapolation [S] to infinite basis set for Z=3 



Species 



Lii- 
Be 

Bi+ 

Ne6+ 
Mg8+ 
Cai6+ 

Ni24+ 



-So 



-Eq; 



-Evjf, 



Iconf 



7.500758 7.49812(6) 7.5008(1 

14.66736 14.65715(4) 14.66707(7 

24.34893 24.3300(2) 24.3486(1 

110.29092 110.2167(2) 110.2902(1 

162.17108 162.073(2) 162.170(1 

469.69474 469.503(7) 469.6935(1 

937.21951 936.936(3) 937.217(9 



Table 2 

Expectation values of radius, r, for one particle numerical Hartree-Fock orbitals 

given in bohr 

Z 

3 

4 

5 

10 

12 

20 

28 



{r)is 


{r)2s 


0.572876 


3.87342 


0.414896 


2.64902 


0.325176 


1.7977 


0.156101 


0.710569 


0.122388 


0.571132 


0.0765909 


0.32221 


0.0515183 


0.207193 



Increasing Z from 3 to 28 has the effect of localizing the single particle orbitals 
radially with corresponding increase in the density closer to nucleus. This can 
be measured by the expectation value of the radius for the orbitals {r)is and 
{r)2s and each change by an order of magnitude (shown in Table |2]). This 
clearly reveals that with increasing density the fixed-node error grows since 
the HF nodes are static and they do not depend on the density (we remind 
the HF node equation (ri — r2){r^ — r^) = 0. 

This is confirmed by the small fixed- node error of the Z = 3 case, less than 
4% of Ecorr, which is much smaller than the error for Be which is larger than 
10% of Ecorr- This is in line with the fact that the Li~^ 2s electrons are highly 
diffuse and Jastrow-type correlations play much more important role overall. 



One way to look at the nodes is to scan the wave function with a "probe" 
electron or pair of electrons for a given position for the rest of the particles. 
One gets a view of the subspace of the 11- dimensional node which corresponds 
to our system. We therefore fixed one of the spin-up electrons at {r, 6, 0} = 
{1.0, f|, ^} and similarly, one of the spin-down electrons at {0.985, |, |}. We 
then scanned the space by the remaining pair of electrons and we draw the 
zero isosurface of the '^2conf, see Fig. [21 The scan shows clearly the effect of 
two configurations which enable the pair of the electrons to "visit" both inside 
and outside region of the two almost concentric spheres which are fused on 
one side creating thus a nodal "opening." This can be contrasted with the HF 
node which is given exactly by two ideal, concentric spheres. During the QMC 
walker evolution process the pair of particles thus cannot get from outside to 
inside-it will be forever locked in one of the (four) equivalent domains. 

These results clearly indicate that the impact of the density on the fixed-node 
errors is very significant. In particular, comparison of Z = 28 with Z = 3 is 
perhaps the most revealing. For Z = 28 the electron density is highly localized 
around the nucleus, and the HF orbitals are very close to noninteracting ones 
since the electron-electron interaction is much smaller than the electron-ion 
contribution. However, the nodal error is significant both in absolute terms 
(283 mHa) and in relative terms (28.6 % of Ecorr)- Interestingly, this system 
is considered "easy" in the sense that the one-particle contributions dominate 
overall. On the other hand, Li~ exhibits nontrivial many-body correlations 
related to the weak anionic bonding which is usually more demanding to 
describe by traditional approaches based on basis sets and excitations. In 
fact, the FNDMC result is surprisingly good even with the HF nodes, which, 
as explained above, have incorrect topology. This seems counterintuitive but 
can be understood once we realize that significantly lower density of Li~ make 
the nodal errors less important. 

The second interesting observation is that the two-configuration wave function 
is very accurate and captures nodes correctly in the range of systems, from the 
weakly bonded anion to the very highly ionized cations. Although it is easy 
to understand that the near-degeneracy effect is crucial it is interesting to 
see that higher excitations appear to have only marginal effect for all studied 
systems. 



4 Conclusions 



We studied the dependence of the fixed-node error on electron density. For this 
purpose we have chosen the system of four electrons in a Coulomb potential 
with varying nuclear charge Z. This system has somewhat unique properties: 
the HF wave function's fixed-node bias is rather large, the HF node is invariant 




Fig. 2. 3D subspace of the 2-configuration nodal surface in real space. The two dots 
at the opening represent the spin-up and -down electrons fixed at slightly different 
radial distances. The tiny dark spot in the middle is the nucleus. The node is found 
by scanning the space with the remaining two electrons located on the top of each 
other and plotting the wave function's zero isosurface. 3 lighting sources are used 
to make the curvature of the surface visible. The semi-transparency enables to see 
"inside" and show that the pair of the scanning electrons can sample both inside 
and outside regions by passing through the opening (i.e., without crossing a node). 
This is not the case for the HF wave function which has the nodal surface always as 
two concentric ideal spheres (one corresponding to spin-up the other to spin-down 
subspaces). 

to varying Z, and for the Be atom it was known that very accurate nodes were 
not too difficuh to construct. Expanding on these insights, our calculations 
have convincingly shown: 



a) the fixed-node bias increases with increasing density and appears to be 
predominantly linear with Z and for the case of four-electron Be-like Coulomb 
system is very closely related to the first-order correction in the MP2 theory; 

b) the 2-configuration wave function produces highly accurate nodes across 



the values of Z we tested, i.e. it appears to be very consistent and basically 
universal for this system. 

This has important implications for the further analysis of the fixed-node 
errors for larger and more complicated systems since the results suggest that 
even small fixed-node errors in the region with high electronic density could 
have much stronger influence on the total correlation energy recovered than 
large fixed-node errors in low electronic density regions. 
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